In [1]:
%pylab inline
import scipy
import scipy.fftpack
import scipy.signal
import scipy.io
load_path = "/home/jbm/capture_commercial.raw"
Fc = 463.4e6
Fs = 2e6
N = -1 # Number of samples to pull
pylab.rcParams['figure.figsize'] = (15.0, 6.0)
In [2]:
x = np.fromfile(load_path, dtype="complex64", count=N)
samps = [ 512000, 128000, 1024000, 16000 ]
#x = sdr.read_samples(2*sdr.sample_rate)
N = 4
figure(figsize=(15,10))
for i,s in enumerate(samps):
subplot(N/2, 2, i+1)
psd(x[:s], Fs=Fs)
title(s)
None
In [3]:
samps = [ 16000, 128000, 512000, 1024000 ]
N = len(samps)
figure(figsize=(15,10))
for i,s in enumerate(samps):
subplot(N/2, 2, i+1)
fs, b = scipy.signal.welch(x[0:s], fs=Fs, nperseg=1024)
fs = scipy.fftpack.fftshift(fs)
b = scipy.fftpack.fftshift(b)
plot(fs, np.log(b))
xlim(-Fs/2, Fs/2)
title(s)
None
In [4]:
fs, b = scipy.signal.welch(x, fs=Fs, nperseg=1024)
fs = scipy.fftpack.fftshift(fs)
b = scipy.fftpack.fftshift(b)
bl = np.log(b)
lpf = scipy.signal.remez(40, [0, 0.001, 0.005, 0.5], [1, 0])
bp = scipy.signal.lfilter(lpf, 1.0, bl)
plot(bp, linewidth = 2)
#ylim(-30, -10)
figure()
plot((bl - bp)[50:])
Out[4]:
Okay, trying to remove the floor with averaging looks like a losing proposition at this point. I could probably better align the two datasets, but let's see what straight-up diffs look like.
In [5]:
plot(diff(np.log(b)))
plot((np.log(b) + 30)/5 )
Out[5]:
In [6]:
hist(diff(np.log(b)), bins=100); None
bl = np.log(b)
dbl = diff(bl)
[ f(dbl) for f in [ np.mean, np.std] ]
Out[6]:
In [7]:
plot(dbl)
axhline(np.mean(dbl) + np.std(dbl)*1.5)
None
This misses a bunch of stuff,because we have a bunch of outliers. Let's focus on the IQR, maybe?
In [8]:
y = sorted(dbl)
[ y[i*len(y)/100] for i in [10, 25, 50, 75, 90] ]
Out[8]:
In [9]:
il = 5*len(y)/100
ih = 95*len(y)/100
yp = y[il:ih]
[ f(yp) for f in [ np.mean, np.std] ]
Out[9]:
In [10]:
plot(dbl)
my_floor = np.mean(yp) + np.std(yp)*3.0
axhline(my_floor)
Out[10]:
In [11]:
figure(figsize=(15,8))
plot(fs, bl)
hits = (np.array(dbl) > my_floor).nonzero()[0]+1
plot(fs[hits], bl[hits], 'o')
xlim(-Fs/2, Fs/2)
Out[11]:
Okay, I like this heuristic. Now to figure out channel widths. Let's ignore the top and bottom 10% of the band for now.
In [12]:
hits = hits[ (hits > len(bl)/10) & (hits < len(bl)*9/10) ]
In [13]:
# Do some janky RLEish work to find contiguous channels.
# There are better numpy-ish ways to do this, but this is legible.
runs = []
last_start = None
gap_allowed = 2 # Allow some fudge factor here
for (h0,h1) in zip(hits[:-1], hits[1:]):
if h1 - h0 < gap_allowed:
if None == last_start:
last_start = h0
else:
if None != last_start:
runs.append( (last_start, h0) )
last_start = None
plot(fs, bl)
for h0,h1 in runs:
axvspan(fs[h0-1], fs[h1+1], color="red", alpha=0.2)
xlim(fs[len(bl)/10], fs[9*len(bl)/10])
[ ((fs[h0]+fs[h1])/2, fs[h1]-fs[h0], fs[h0], fs[h1]) for (h0,h1) in runs ]
Out[13]:
